Language/OS - Multiplatform Resource Library/LANGUAGE OS.iso/a_utils/expanded.lha/expanded/src/Space.C

C/C++ Source or Header | 1992-03-19 | 71KB | 2,597 lines

//
// Linear-Affine-Projective Geometry Package
//
// Space.C
//
// $Header$
//
// William J.R. Longabaugh
// University of Washington
//
// Implementation of the linear-affine-projective geometry
// package described in William J.R. Longabaugh, "An Expanded
// System for Coordinate-Free Geometric Programming", Master's
// thesis, University of Washington, 1992.
//
// Copyright (c) 1992, William J.R. Longabaugh
//   Copying, use and development for non-commercial purposes permitted.
//   All rights for commercial use reserved.
//   This software is unsupported and without warranty; it is
//   provided "as is".
//
// ***********************************************************************

#include <string.h>
#include "Lap.h"

// ***********************************************************************
//
// Class for exclusive use of Space class to hold onto data for a
// space.  What is really wanted is a simple structure; I am not
// really interested in data hiding here.  However, under g++ 1.37.1,
// this apparently MUST be a class, not a struct, so that virtual
// ptrs of class members are initialized.  Also, the fact that class
// members will be automatically initialized using null constructors
// is desirable.  Thus, this is a class, but note that all members
// are public, allowing unlimited access.

// Constants to organize space set and standard affine map arrays:

const int SP_SET_SIZE = 6;

const int NUM_SLOT = 3;
const int AFFSLOT = 0;
const int LINSLOT = 1;
const int PCSLOT = 2;

class SpaceInfo {

  public:

    SpaceType type;            // What type of space it is
    SRel thisSpace;            // Position of this space in 6-set
    int dim;             // The dimension of the space
    char name[MAX_NAMESIZE];     // Space name
    Basis stdBasis;         // Predefined Basis
    int matsize;            // Size of matrix reps
    GeObType native;        // Native GeOb type

    int cpSize;            // Number of spaces, if a CPSpace
    SpaceList cpSpaces;        // Constituent spaces (if CPSize > 1)
    IntList cpSlots;        // Starting slots for components

    Space spaceSet[SP_SET_SIZE];     // Other spaces in space set
    VSpace dual;             // Dual space

    ProjectiveMap stdPMap;        // Standard Projective Map
        AffineMap stdAMap[NUM_SLOT];     // Standard Affine Maps
    LinearMap stdLMap;        // Standard Linear Map

    Boolean     isEuclidean;        // True if space is Euclidean
    MLM      innerProd;        // Standard Inner Product
    MLM     crossProd;        // Standard Cross Product

    GeObList atInfinity;        // Set of points at infinity.
    Matrix     hoistTangent;         // Gets tan. vec. to affine rep.

    ASubSet stdAffineSet;         // Std. affine subset of space
    SubSet fullset;            // SubSet for whole space
    PSubSet fullproj;        // Peoj. SubSet for whole space

// Constructors/Destructors:

        SpaceInfo(void);          // Null constructor
        ~SpaceInfo(void);         // Destructor
        SpaceInfo(char* namein, int d); // Used to build error info block
};

// ***********************************************************************
//
// Declare the existence of the error info block:

static SpaceInfo errinfo;

// ***********************************************************************
// ***********************************************************************
//
// SpaceInfo Class
//
// ***********************************************************************
// ***********************************************************************
//
// Null constructor.  Nothing needs to be done.  Fortunately, the
// null constructors of object class members are also called automatically,
// so they are all initialized to null objects (which is important).
// This is too large to inline in g++.
//

SpaceInfo::SpaceInfo(void) {}

// ***********************************************************************
//
// Destructor.  Don't actually need this, but it suppresses a g++ 1.37
// compiler warning.  It is also too large to inline.

SpaceInfo::~SpaceInfo(void) {}

// ***********************************************************************
//
// Special constructor used to build an "error info block".  Uninitialized
// spaces point to that block.

SpaceInfo::SpaceInfo(char* namein, int d)
{
  type = NULL_SPACE;
  dim = d;
  thisSpace= NO_RELATION;
  matsize = 0;
  native = NULL_GEOB;
  cpSize = 0;
  isEuclidean = FALSE;

// Local copy of the debug name:

  strncpy(name, namein, MAX_NAMESIZE - 1);
  name[MAX_NAMESIZE - 1] = '\0';
}

// ***********************************************************************
// ***********************************************************************
//
// Space Class
//
// ***********************************************************************
// ***********************************************************************
//
// Private/protected member functions
//
// ***********************************************************************
//
// NOTE: The following few functions are PERFECT candidates for inlining.
// However, because they reference items in a SpaceInfo, they cannot be
// defined in the Space class declaration.  (The SpaceInfo class must
// be declared following the List classes, which must in turn be declared
// after the Space class.)  While the "inline" keyword could be used
// to declare the function inline when they are defined, the CFRONT 1.2
// compiler appears to have a bug that occasionally shows up when this
// is done.  Because of this problem, it was decided that the only inline
// functions in this system would be those defined in the class
// declaration.
//
// ***********************************************************************
//
// Note it is the responsibility of caller to check validity of bounds:

int Space::StartSlot(int n) {return (info->cpSlots[n]);}

// ***********************************************************************

int Space::MatrixSize(void) {return (info->matsize);}

// ***********************************************************************

Matrix& Space::HoistTangent(void) {return (info->hoistTangent);}

// ***********************************************************************
//
// User is responsible for updating space set independently!

inline void Space::SetDual(Space& s) {info->dual = s;}

// ***********************************************************************

void Space::SetStdPMap(Map& m) {info->stdPMap = m;}

// ***********************************************************************

void Space::SetStdLMap(Map& m) {info->stdLMap = m;}

// ***********************************************************************

void Space::SetStdAMap(SRel s, Map& m)
{
  switch (s) {
    case LINEARIZATION:
      info->stdAMap[LINSLOT] = m;
      break;
    case PROJECT_COMP:
      info->stdAMap[PCSLOT] = m;
      break;
    case AFFINE:
      info->stdAMap[AFFSLOT] = m;
      break;
    case TANGENT:
    case LIN_DUAL:
    case TANG_DUAL:
    case NO_RELATION:
    case SAME_SPACE:
    default:
      errh.ErrorExit("void Space::SetStdAMap(SRel, Map&)",
                     "Unexpected Relation",
              ErrType("Space relation specified =", s, SRELTYPES),
              *this, m);
      break;
  }

  info->stdAffineSet = m.Domain();

  if (info->thisSpace == PROJECT_COMP) {
    info->atInfinity = m.Domain().AtInfinity();
  }
}

// ***********************************************************************

void Space::SetSpaceSet(SRel r, Space& sp)
{
  if (r >= (int)NO_RELATION) {
    errh.ErrorExit("void Space::SetSpaceSet(SRel, Space&)",
                   "Invalid relation specified",
            ErrType("Space relation specified =", r, SRELTYPES),
            *this, sp);
  }
  info->spaceSet[(int) r] = sp;
}

// ***********************************************************************
//
// This private method broadcasts the existence of a newly created space
// to all other spaces in the set.  It assumes the set of the new space
// has been initialized to null spaces (except for its own slot).  The
// SRel argument refers to the relation of this space.
//

void Space::Broadcast(SRel s, Space& ins)
{
// The space we were handed gives all the spaces that are related to this
// space.  We need to fill this information in, AND broadcast the existence
// of this new space to the other spaces in the set:

  for (int i = 0; i < SP_SET_SIZE; i++) {
    if ((SRel)i == s) {
      ins.SetSpaceSet(s, *this);
    } else if (ins.HasSpace((SRel)i)) {
      Space temp = ins.GetSpace((SRel)i);
      info->spaceSet[i] = temp;
      temp.SetSpaceSet(s, *this);
    }
  }
}

// ***********************************************************************
//
// This function merges two independent space sets into a single set.
//

void Space::MergeSets(Space& set2)
{

// Make this space's set up to date by merging the two lists.  Then
// for everybody on this new list (except this space) tell them about
// all the spaces on the merged list:

  int i;
  Space hold;

  for (i = 0; i < SP_SET_SIZE; i++) {
    if (set2.HasSpace((SRel)i)) {
      info->spaceSet[i] = set2.GetSpace((SRel)i);
    }
  }

  for (i = 0; i < SP_SET_SIZE; i++) {
    hold = info->spaceSet[i];
    if ((hold != *this) && (hold.Holds() != NULL_SPACE)) {
      for (int j = 0; j < SP_SET_SIZE; j++) {
        if (!hold.HasSpace((SRel)j)) {
      hold.SetSpaceSet((SRel)j, info->spaceSet[j]);
        }
      }
    }
  }
}

// ***********************************************************************
//
// This function copies only a maximum number of characters into
// the debug name buffer.
//

void Space::CopyDebugName(char* buf)
{
  strncpy(info->name, buf, MAX_NAMESIZE - 1);
  info->name[MAX_NAMESIZE - 1] = '\0';
  return;
}

// ***********************************************************************
//
// This function builds a tagged debug name.  The user is responsible
// for calling DestroyTagName() when finished.
//

char* Space::BuildTagName(char* tag, Space& s)
{
  int len = strlen(tag);
  char* buf = new char[MAX_NAMESIZE + len];

  strncpy(buf, tag, len);
  s.Name(buf + len);
  return (buf);
}

// ***********************************************************************
//
// Destroys a tagged debug name:
//

void Space::DestroyTagName(char* buf) {delete buf;}

// ***********************************************************************
//
// Evaluate the linkage between two given spaces.  They are unlinked,
// linked to each other, or erroneously linked to other spaces:
//

Boolean Space::HaveLink(Space& s1, SRel r1,
            Space& s2, SRel r2, Boolean* have_error)
{
  Boolean linked = FALSE;
  Boolean s2linked = FALSE;
  Space s1c;
  Space s2c;

  if (s1.HasSpace(r2)) {
    s1c = s1.GetSpace(r2);
    linked = TRUE;
  }

  if (s2.HasSpace(r1)) {
    s2c = s2.GetSpace(r1);
    s2linked = TRUE;
  }

  if ((linked && (s1c != s2)) || (s2linked && (s2c != s1))) {
    *have_error = TRUE;
  } else {
    *have_error = FALSE;
  }

  return linked;
}

// ***********************************************************************
//
// The following group of functions are used to link spaces together
// with standard maps.  A variety of functions are needed to handle
// the many different cases of how the spaces may be already linked.
//
// ***********************************************************************
//
// Links the Tangent space to the Linearization.  The LA link must
// already exist (or the tangent space would not exist), so the
// map between them exists.
//

void Space::LinkLT(void)
{
  Space T = this->GetSpace(TANGENT);
  Space L = this->GetSpace(LINEARIZATION);
  Space A = this->GetSpace(AFFINE);

  LinearMap m = (A.AffineMapTo(LINEARIZATION)).AssocLinear();
  T.SetStdLMap(m);
  L.SetStdLMap(m.Inv());
  return;
}

// ***********************************************************************
//
// Links the Affine space to the Projective Completion, where there
// are no preexisting links:
//

void Space::LinkAP(AffineMap& m)
{
  Space A = this->GetSpace(AFFINE);
  Space P = this->GetSpace(PROJECT_COMP);

  A.SetStdAMap(PROJECT_COMP, m);
  P.SetStdAMap(AFFINE, m.Inv());
  return;
}

// ***********************************************************************
//
// Links the Linearization space to Affine space, where there
// are no preexisting links:
//

void Space::LinkLA(AffineMap& m)
{
  Space L = this->GetSpace(LINEARIZATION);
  Space A = this->GetSpace(AFFINE);
  AffineMap mi = m.Inv();

  L.SetStdAMap(AFFINE, m);
  A.SetStdAMap(LINEARIZATION, mi);

  if (this->HasSpace(TANGENT)) {
    Space T = this->GetSpace(TANGENT);
    LinearMap mia = mi.AssocLinear();
    T.SetStdLMap(mia);
    L.SetStdLMap(mia.Inv());
  }
  return;
}

// ***********************************************************************
//
// Links the Linearization space to the Projective Completion, where
// there are no preexisting links:
//

void Space::LinkLP(ProjectiveMap& m)
{
  Space L = this->GetSpace(LINEARIZATION);
  Space P = this->GetSpace(PROJECT_COMP);

  L.SetStdPMap(m);
  P.SetStdPMap(m.Inv());
  return;
}

// ***********************************************************************
//
// Links the Linearization space to Affine space, where there
// is a link from Lin to PC:
//

void Space::LinkLAHaveLP(AffineMap& m)
{
  Space L = this->GetSpace(LINEARIZATION);
  Space A = this->GetSpace(AFFINE);
  Space P = this->GetSpace(PROJECT_COMP);
  AffineMap mi = m.Inv();

  L.SetStdAMap(AFFINE, m);
  A.SetStdAMap(LINEARIZATION, mi);

  if (this->HasSpace(TANGENT)) {
    Space T = this->GetSpace(TANGENT);
    LinearMap mia = mi.AssocLinear();
    T.SetStdLMap(mia);
    L.SetStdLMap(mia.Inv());
  }

// Since a standard projective map exists, the AP map must be derived:

  AffineMap ap(L.ProjectiveMapTo(PROJECT_COMP), m);

  A.SetStdAMap(PROJECT_COMP, ap);
  P.SetStdAMap(AFFINE, ap.Inv());

// Finally compose the affine maps to get the LP affine map:

  AffineMap lp = ap.Compose(m);

  L.SetStdAMap(PROJECT_COMP, lp);
  P.SetStdAMap(LINEARIZATION, lp.Inv());
  return;
}

// ***********************************************************************
//
// Links the Affine space to the Projective Completion, where there
// is a link from Lin to PC:
//

void Space::LinkAPHaveLP(AffineMap& m)
{
  Space L = this->GetSpace(LINEARIZATION);
  Space A = this->GetSpace(AFFINE);
  Space P = this->GetSpace(PROJECT_COMP);

  A.SetStdAMap(PROJECT_COMP, m);
  P.SetStdAMap(AFFINE, m.Inv());

// Standard projective PL map exists, so LA map must be derived:

  AffineMap la(m, L.ProjectiveMapTo(PROJECT_COMP));

  A.SetStdAMap(LINEARIZATION, la.Inv());
  L.SetStdAMap(AFFINE, la);

// Compose the affine maps to get the LP affine map:

  AffineMap lp = m.Compose(L.AffineMapTo(AFFINE));

  L.SetStdAMap(PROJECT_COMP, lp);
  P.SetStdAMap(LINEARIZATION, lp.Inv());

// Fill in standard linear map, if tangent space exists:

  if (this->HasSpace(TANGENT)) {
    Space T = this->GetSpace(TANGENT);
    LinearMap ala = la.Inv().AssocLinear();
    T.SetStdLMap(ala);
    L.SetStdLMap(ala.Inv());
  }
  return;
}

// ***********************************************************************
//
// Links the Linearization space to Affine space, where there
// is a link from Aff to PC:
//

void Space::LinkLAHaveAP(AffineMap& m)
{
  Space L = this->GetSpace(LINEARIZATION);
  Space A = this->GetSpace(AFFINE);
  Space P = this->GetSpace(PROJECT_COMP);
  AffineMap mi = m.Inv();

  L.SetStdAMap(AFFINE, m);
  A.SetStdAMap(LINEARIZATION, mi);

  if (this->HasSpace(TANGENT)) {
    Space T = this->GetSpace(TANGENT);
    LinearMap mia = mi.AssocLinear();
    T.SetStdLMap(mia);
    L.SetStdLMap(mia.Inv());
  }

// Compose the affine maps to get the LP affine map:

  AffineMap lp = A.AffineMapTo(PROJECT_COMP).Compose(m);

  L.SetStdAMap(PROJECT_COMP, lp);
  P.SetStdAMap(LINEARIZATION, lp.Inv());

// Now derive a consistent projective map for lin-pc link:

  ProjectiveMap plp(L, lp, P);

  L.SetStdPMap(plp);
  P.SetStdPMap(plp.Inv());
  return;
}

// ***********************************************************************
//
// Links the Affine space to the Projective Completion, where there
// is a link from Lin to Aff:
//

void Space::LinkAPHaveLA(AffineMap& m)
{
  Space L = this->GetSpace(LINEARIZATION);
  Space A = this->GetSpace(AFFINE);
  Space P = this->GetSpace(PROJECT_COMP);

  A.SetStdAMap(PROJECT_COMP, m);
  P.SetStdAMap(AFFINE, m.Inv());

// Compose the affine maps to get the LP affine map:

  AffineMap lp = m.Compose(L.AffineMapTo(AFFINE));

  L.SetStdAMap(PROJECT_COMP, lp);
  P.SetStdAMap(LINEARIZATION, lp.Inv());

// Now derive a consistent projective map for lin-pc link:

  ProjectiveMap plp(L, lp, P);

  L.SetStdPMap(plp);
  P.SetStdPMap(plp.Inv());
  return;
}

// ***********************************************************************
//
// Links the Linearization space to the Projective Completion, where
// is a link from Lin to Aff:
//

void Space::LinkLPHaveLA(ProjectiveMap& m)
{
  Space L = this->GetSpace(LINEARIZATION);
  Space A = this->GetSpace(AFFINE);
  Space P = this->GetSpace(PROJECT_COMP);

  L.SetStdPMap(m);
  P.SetStdPMap(m.Inv());

// Since a standard projective map now exists, the AP map must be derived:

  AffineMap ap(m, L.AffineMapTo(AFFINE));

  A.SetStdAMap(PROJECT_COMP, ap);
  P.SetStdAMap(AFFINE, ap.Inv());

// Finally compose the affine maps to get the LP affine map:

  AffineMap lp = ap.Compose(L.AffineMapTo(AFFINE));

  L.SetStdAMap(PROJECT_COMP, lp);
  P.SetStdAMap(LINEARIZATION, lp.Inv());
  return;
}

// ***********************************************************************
//
// Links the Linearization space to the Projective Completion, where
// is a link from Aff to PC:
//

void Space::LinkLPHaveAP(ProjectiveMap& m)
{
  Space L = this->GetSpace(LINEARIZATION);
  Space A = this->GetSpace(AFFINE);
  Space P = this->GetSpace(PROJECT_COMP);
  ProjectiveMap mi = m.Inv();

  L.SetStdPMap(m);
  P.SetStdPMap(mi);

// Standard affine AP map exists, so LA map must be derived:

  AffineMap ap = A.AffineMapTo(PROJECT_COMP);
  AffineMap la(ap, m);

  A.SetStdAMap(LINEARIZATION, la.Inv());
  L.SetStdAMap(AFFINE, la);

// Compose the affine maps to get the LP affine map:

  AffineMap lp = ap.Compose(L.AffineMapTo(AFFINE));

  L.SetStdAMap(PROJECT_COMP, lp);
  P.SetStdAMap(LINEARIZATION, lp.Inv());

// Fill in standard linear map, if tangent space exists:

  if (this->HasSpace(TANGENT)) {
    Space T = this->GetSpace(TANGENT);
    LinearMap ala = la.Inv().AssocLinear();
    T.SetStdLMap(ala);
    L.SetStdLMap(ala.Inv());
  }
  return;
}

// ***********************************************************************
//
// Returns a REFERENCE to the standard affine map to the specified
// space in the six-space set.  Only used by the implementing classes.
//

AffineMap& Space::AffineMapToRef(SRel s)
{
  static char* sig = "AffineMap Space::AffineMapToRef(SRel)";
  Boolean errval = FALSE;

// Need to get the space first to insure it is created.  Note if
// we are trying to get from linearization to projective with affine
// map, the intervening affine space must be created also.  LIN to PC
// links without an affine space are purely projective; standard
// affine map is only created when affine space is plugged in.

  if (((info->thisSpace == LINEARIZATION) && (s == PROJECT_COMP)) ||
      ((info->thisSpace == PROJECT_COMP) && (s == LINEARIZATION))) {
    this->GetSpace(AFFINE);
  }

  this->GetSpace(s);

  if (info->thisSpace == LINEARIZATION) {
    if (s == AFFINE) {
      return info->stdAMap[AFFSLOT];
    } else if (s == PROJECT_COMP) {
      return info->stdAMap[PCSLOT];
    } else {
      errval = TRUE;
    }
  } else if (info->thisSpace == AFFINE) {
    if (s == LINEARIZATION) {
      return info->stdAMap[LINSLOT];
    } else if (s == PROJECT_COMP) {
      return info->stdAMap[PCSLOT];
    } else {
      errval = TRUE;
    }
  } else if (info->thisSpace == PROJECT_COMP) {
    if (s == LINEARIZATION) {
      return info->stdAMap[LINSLOT];
    } else if (s == AFFINE) {
      return info->stdAMap[AFFSLOT];
    } else {
      errval = TRUE;
    }
  } else {
    errh.ErrorExit(sig, "No standard affine map from this space", *this);
  }

  if (errval) {
    errh.ErrorExit(sig, "No standard affine map to specified space",
           ErrType("Space relation specified =", s, SRELTYPES),
           *this);
  }
}

// ***********************************************************************
//
// Returns a REFERENCE to the standard projective map to the specified
// space in the six-space set.  Only used by the implementing classes.
//

ProjectiveMap& Space::ProjectiveMapToRef(SRel s)
{

// Need to get the space first to insure it is created!

  this->GetSpace(s);

  if (((info->thisSpace == LINEARIZATION) && (s == PROJECT_COMP)) ||
      (info->thisSpace == PROJECT_COMP) && (s == LINEARIZATION)) {
    return info->stdPMap;
  } else {
    errh.ErrorExit("ProjectiveMap Space::ProjectiveMapToRef(SRel)",
                   "No standard projective map to specified space",
            ErrType("Space relation specified =", s, SRELTYPES),
           *this);
  }
}

// ***********************************************************************
//
// Returns a REFERENCE to the standard linear map to the specified
// space in the six-space set.  Only used by the implementing classes.
//

LinearMap& Space::LinearMapToRef(SRel s)
{
// Need to get the space first to insure it is created!

  this->GetSpace(s);

  if (((info->thisSpace == LINEARIZATION) && (s == TANGENT)) ||
      (info->thisSpace == TANGENT) && (s == LINEARIZATION)) {
    return info->stdLMap;
  } else {
    errh.ErrorExit("LinearMap Space::LinearMapToRef(SRel)",
                   "No standard linear map to specified space",
            ErrType("Space relation specified =", s, SRELTYPES),
           *this);
  }
}

// ***********************************************************************
//
// Data dump for short form.
//

void Space::data_out(ostream &c, int indent, char* name)
{
  char *ibloc = new char[indent + 1];
  for (int i = 0; i < indent; i++) {
    *(ibloc + i) = ' ';
  }
  *(ibloc + indent) = '\0';

  c << ibloc << ast;
  c << ibloc << name;
  c << ibloc << "Type of space: ";
  SpaceTypeOut(c, type);
  c << "\n";
  c << ibloc << "Currently holds: ";
  SpaceTypeOut(c, holds);
  c << "\n";
  if (holds != NULL_SPACE) {
    c << ibloc << "Space relationship: ";
    SRelOut(c, info->thisSpace);
    c << "\n" << ibloc << "Name of space: " << info->name << "\n";
    c << ibloc << "Dimension: " << info->dim << "\n";
    c << ibloc << "Size of matrix representations: " << info->matsize << "\n";
    c << ibloc << "Native type of object in space: ";
    GeObTypeOut(c, info->native);
    c << "\n";
    c << ibloc << "Memory location: " << (long)info << "\n";
    c << ibloc << "Number of constituent spaces: " << info->cpSize << "\n";

    if (info->cpSize > 1) {
      c << ibloc << "Constituent spaces:\n";
      for (i = 0; i < info->cpSize; i++) {
        (info->cpSpaces[i]).debug_out(c, indent + ERR_IND);
      }
    }
  }
  c << ibloc << ast;

  delete ibloc;
  return;
}

// ***********************************************************************
//
// This dumps ALL the data in the space.  It should not be
// be called by the debug output functions of any class that has a
// member in space, or infinite recursion will result.
//

void Space::heavy_data_out(ostream &c, int indent, char* name)
{
  int i;

  char *ibloc = new char[indent + 1];
  for (i = 0; i < indent; i++) {
    *(ibloc + i) = ' ';
  }
  *(ibloc + indent) = '\0';

  c << ibloc << ast;
  c << ibloc << name;
  c << ibloc << "Type of space: ";
  SpaceTypeOut(c, type);
  c << "\n";
  c << ibloc << "Currently holds: ";
  SpaceTypeOut(c, holds);
  c << "\n";
  if (holds != NULL_SPACE) {
    c << ibloc << "Space relationship: ";
    SRelOut(c, info->thisSpace);
    c << "\n" << ibloc << "Name of space: " << info->name << "\n";
    c << ibloc << "Dimension: " << info->dim << "\n";
    c << ibloc << "Memory location: " << (long)info << "\n";
    c << ibloc << "Standard basis:\n";
    info->stdBasis.debug_out(c, indent + ERR_IND);
    c << ibloc << "Size of matrix representations: " << info->matsize << "\n";
    c << ibloc << "Native type of object in space: ";
    GeObTypeOut(c, info->native);
    c << "\n";
    c << ibloc << "Number of constituent spaces: " << info->cpSize << "\n";

    if (info->cpSize > 1) {
      c << ibloc << "Constituent spaces:\n";
      for (i = 0; i < info->cpSize; i++) {
        (info->cpSpaces[i]).debug_out(c, indent + ERR_IND);
      }

      c << ibloc << "Starting locations for each constituent:\n";
      for (i = 0; i < info->cpSize; i++) {
        c << ibloc << info->cpSlots[i];
        if (i % 4 == 3) {c << "\n";}
      }
    }

    c << ibloc << "Other spaces in space set:\n";
    for (i = 0; i < 6; i++) {
      if ((info->spaceSet[i]).Holds() != NULL_SPACE) {
        c << "\n" << ibloc;
        SRelOut(c, i);
        c << ":\n";
        (info->spaceSet[i]).debug_out(c, indent + ERR_IND);
      }
    }

    if ((info->dual).Holds() != NULL_SPACE) {
      c << ibloc << "Dual Space:\n";
      info->dual.debug_out(c, indent + ERR_IND);
    }

    c << ibloc << "Standard projective map:\n";
    info->stdPMap.debug_out(c, indent + ERR_IND);
    c << ibloc << "Standard affine maps:\n";
    for (i = 0; i < NUM_SLOT; i++) {
      info->stdAMap[i].debug_out(c, indent + ERR_IND);
    }
    c << ibloc << "Standard linear map:\n";
    info->stdLMap.debug_out(c, indent + ERR_IND);

    if (info->isEuclidean == TRUE) {
      c << ibloc << "\nIs a Euclidean space\n";
      c << ibloc << "Standard inner product:\n";
      info->innerProd.debug_out(c, indent + ERR_IND);
      c << ibloc << "Standard cross product:\n";
      info->crossProd.debug_out(c, indent + ERR_IND);
    } else {
      c << ibloc << "\nIs not a Euclidean space\n";
    }

    c << ibloc << "Set of points at infinity:\n";
    info->atInfinity.debug_out(c, indent + ERR_IND);

    c << ibloc << "Tangent vector hoisting matrix:\n";
    info->hoistTangent.debug_out(c, indent + ERR_IND);

    c << ibloc << "Standard affine subset:\n";
    info->stdAffineSet.debug_out(c, indent + ERR_IND);

    c << ibloc << "Full space subset:\n";
    info->fullset.debug_out(c, indent + ERR_IND);

    c << ibloc << "Full space projective subset:\n";
    info->fullproj.debug_out(c, indent + ERR_IND);
  }
  c << ibloc << ast;

  delete ibloc;
  return;
}
// ***********************************************************************


// ***********************************************************************
//
// Public member functions
//
// ***********************************************************************

Space::Space(void) {type = ANY_SPACE; holds = NULL_SPACE; info = &errinfo;}

// ***********************************************************************
//
// Do memberwise initialization.  Since no class members need to be
// copied, this constructor should not be necessary.  However,
// it seems to be necessary to define this to suppress a wierd compiler
// error.
//

Space::Space(Space& v) {type = ANY_SPACE; holds = v.holds; info = v.info;}

// ***********************************************************************
//
// Do memberwise assignment.  (Needed to catch g++ 1.37 bug).
//

Space& Space::operator=(Space &s)
{
//
// This weird checking is required to bypass the apparent inheritance of
// assignment under g++ 1.37:
//
  if ((type != ANY_SPACE) &&
      ((type != s.holds) && (s.holds != NULL_SPACE))) {
    errh.ErrorExit("Space& Space::operator=(Space &)",
                   "Illegal assignment attempted", *this, s);
  }
  holds = s.holds;
  info = s.info;
  return (*this);
}

// ***********************************************************************
//
// More perfect candidates for inlining, but see the above remark:
//
// ***********************************************************************

int Space::Dim(void) {return (info->dim);}

// ***********************************************************************

char* Space::Name(char *buf) {return (strcpy(buf, info->name));}

// ***********************************************************************

Basis Space::StdBasis(void) {return (info->stdBasis);}

// ***********************************************************************

Simplex Space::StdSimplex(void) {return (Simplex(info->stdBasis));}

// ***********************************************************************

int Space::CPSpaceSize(void) {return (info->cpSize);}

// ***********************************************************************

SRel Space::ThisSpaceIs(void) {return (info->thisSpace);}

// ***********************************************************************

SubSet Space::FullSet(void) {return info->fullset;}

// ***********************************************************************

GeObType Space::NativeType(void) {return (info->native);}

// ***********************************************************************
//
// The following three functions are not inlined to avoid bitwise copying
// under CFRONT 1.2:
//

AffineMap Space::AffineMapTo(SRel s) {return (this->AffineMapToRef(s));}

// ***********************************************************************

ProjectiveMap Space::ProjectiveMapTo(SRel s)
  {return (this->ProjectiveMapToRef(s));}

// ***********************************************************************

LinearMap Space::LinearMapTo(SRel s) {return(this->LinearMapToRef(s));}

// ***********************************************************************

Boolean Space::IsEuclidean(void)
{
  if ((holds != VEC_SPACE) && (holds != AFF_SPACE)) {
    errh.ErrorExit("Boolean Space::IsEuclidean(void)",
           "Space is not vector or affine", *this);
  }
  return (info->isEuclidean);
}

// ***********************************************************************

GeObList Space::PtsAtInfinity(void)
{
  if (holds != PROJ_SPACE) {
    errh.ErrorExit("GeObList Space::PtsAtInfinity(void)",
           "Space is not projective", *this);
  }
  return (info->atInfinity);
}

// ***********************************************************************

PSubSet Space::FullProjSet(void)
{
  if ((holds != VEC_SPACE) && (holds != PROJ_SPACE)) {
    errh.ErrorExit("SubSet Space::FullProjSet(void)",
           "Space is not vector or projective", *this);
  }
  return (info->fullproj);
}

// ***********************************************************************
//
// This is a "short form" debug output for a Space.  It is the form
// that will be called by all other objects that contain a space, and
// will not cause infinite recursion of printouts.  Use the extended
// form to get a full dump of the space data.
//

void Space::debug_out(ostream &c, int indent)
{
  static char *name = "Space object\n";
  this->data_out(c, indent, name);
  return;
}

// ***********************************************************************
//
// This debug output dumps ALL the data in the space.  It should not be
// be called by the debug output functions of any class that has a
// member in space, or infinite recursion will result.
//

void Space::heavy_debug_out(ostream &c, int indent)
{
  static char *name = "Full data dump for a space\n";
  this->heavy_data_out(c, indent, name);
  return;
}


// ***********************************************************************
//
// Given a space, this function returns the SRel indicating how the
// argument is related to this space:
//

SRel Space::SpaceRelation(Space &s)
{
  if (s == *this) {
    return (SAME_SPACE);
  } else {
    for (int i = 0; i < SP_SET_SIZE; i++) {
      if (info->spaceSet[i] == s) {
    return (SRel)i;
      }
    }
    return (NO_RELATION);
  }
}

// ***********************************************************************
//
// Check to see if the related space has been specified, without
// creating it if it has not:
//

Boolean Space::HasSpace(SRel s)
{
  if (s >= (int)NO_RELATION) {
    errh.ErrorExit("Boolean Space::HasSpace(SRel)",
                   "Invalid relation specified",
            ErrType("Space relation specified =", s, SRELTYPES));
  }

// Look in the slot.  If the space we find is not a null space, it has
// been specified, so return TRUE.

  return ((info->spaceSet[s]).Holds() != NULL_SPACE);
}

// ***********************************************************************
//
// Get the space in the specified slot.  If no space of the specified
// relation exists, this routine creates it.
//

Space Space::GetSpace(SRel s)
{
  static char* sig = "Space Space::GetSpace(SRel)";
  Space retval;
  Space holdsp;
  Boolean eflag;
  char* tag;

  if (s >= (int)NO_RELATION) {
    errh.ErrorExit(sig, "Invalid relation specified",
           ErrType("Space relation specified =", s, SRELTYPES));
  }

// Look in the slot.  If the space we find is a null space, it has not
// been created, so build it.

  if ((info->spaceSet[s]).Holds() != NULL_SPACE) {
    return (info->spaceSet[s]);
  } else {
    switch (s) {
      case TANGENT:
      case LIN_DUAL:
      case TANG_DUAL:
    retval = VSpace(s, *this);
        break;
      case LINEARIZATION:
    if (this->HasSpace(AFFINE)) {
      holdsp = this->GetSpace(AFFINE);
      eflag = holdsp.IsEuclidean();
    } else {
      holdsp = this->GetSpace(PROJECT_COMP);
      eflag = TRUE;
    }
    tag = this->BuildTagName("Linearization for ", holdsp);
    retval = VSpace(tag, eflag, holdsp);
    this->DestroyTagName(tag);
        break;
      case PROJECT_COMP:
    if (this->HasSpace(AFFINE)) {
      holdsp = this->GetSpace(AFFINE);
    } else {
      holdsp = this->GetSpace(LINEARIZATION);
    }
    tag = this->BuildTagName("Proj. Space for ", holdsp);
    retval = PSpace(tag, holdsp);
    this->DestroyTagName(tag);
        break;
      case AFFINE:
    if (this->HasSpace(LINEARIZATION)) {
      holdsp = this->GetSpace(LINEARIZATION);
      eflag = holdsp.IsEuclidean();
    } else {
      holdsp = this->GetSpace(PROJECT_COMP);
      eflag = TRUE;
    }
    tag = this->BuildTagName("Affine Space for ", holdsp);
    retval = ASpace(tag, eflag, holdsp);
    this->DestroyTagName(tag);
        break;
      case NO_RELATION:
      case SAME_SPACE:
      default:
        errh.ErrorExit(sig, "Invalid relation specified",
               ErrType("Space relation specified =", s, SRELTYPES));
        break;
      }
  }
  return (retval);
}

// ***********************************************************************
//
// Used to stitch spaces together manually.
//

Space Space::SetSpace(Space& s, Map& m)
{
  static char* sig = "Space Space::SetSpace(Space&, Map&)";
  Space rv;

// Dispatch to different routines depending on what type of space
// this is:

  switch (holds) {
    case VEC_SPACE:
      rv = VSpace(*this).SetSpace(s, m);
      break;
    case AFF_SPACE:
      rv = ASpace(*this).SetSpace(s, m);
      break;
    case PROJ_SPACE:
      rv = PSpace(*this).SetSpace(s, m);
      break;
    case NULL_SPACE:
    case ANY_SPACE:
    default:
      errh.ErrorExit(sig, "Illegal space type",
                 ErrType("Type specified =", holds, SPACETYPES), *this);
      break;
  }
  return (rv);
}

// ***********************************************************************
//
// Returns the dual space
//

VSpace Space::Dual(void)
{
  SRel rel;

// Again, lazy space creation is an issue.

  if ((info->dual).Holds() != NULL_SPACE) {
    return (info->dual);
  } else if (info->thisSpace == LINEARIZATION) {
    rel = LIN_DUAL;
  } else if (info->thisSpace == TANGENT) {
    rel = TANG_DUAL;
  } else {
    errh.ErrorExit("VSpace Space::Dual(void)",
                   "Space cannot have a dual", *this);
  }
  info->dual = VSpace(rel, *this);
  return (info->dual);
}

// ***********************************************************************
//
// Inner products are created lazily:
//

MLM Space::InnerProduct(void)
{
  static char* sig = "MLM Space::InnerProduct(void)";

  if (holds != VEC_SPACE) {
    errh.ErrorExit(sig, "Must be a vector space", *this);
  }
  if (!info->isEuclidean) {
    errh.ErrorExit(sig, "Must be a Euclidean space", *this);
  }

  if (info->innerProd.Holds() == NULL_MULTI) {
    info->innerProd = MLM(INNER_PROD, *this);
  }
  return (info->innerProd);
}

// ***********************************************************************
//
// Cross products are created lazily:
//

MLM Space::CrossProduct(void)
{
  static char* sig = "MLM Space::CrossProduct(void)";

  if (holds != VEC_SPACE) {
    errh.ErrorExit(sig, "Must be a vector space", *this);
  }
  if (!info->isEuclidean) {
    errh.ErrorExit(sig, "Must be a Euclidean space", *this);
  }

  if (info->crossProd.Holds() == NULL_MULTI) {
    info->crossProd = MLM(CROSS_PROD, *this);
  }
  return (info->crossProd);
}

// ***********************************************************************
//
// Only Spaces that have been linked to an affine space have a standard
// affine subset.
//

ASubSet Space::StdAffineSubset(void)
{
  ASubSet retval;

  static char* sig = "ASubSet Space::StdAffineSubset(void)";
  switch (info->thisSpace) {
    case LINEARIZATION:
    case PROJECT_COMP:
      this->GetSpace(AFFINE); // make sure space exists first
      retval = info->stdAffineSet;
      break;
    case AFFINE:
      retval = info->stdAffineSet;
      break;
    case TANGENT:
    case LIN_DUAL:
    case TANG_DUAL:
      errh.ErrorExit(sig,
             "Space does not have a standard affine subset", *this);
      break;
    case NO_RELATION:
    case SAME_SPACE:
    default:
      errh.ErrorExit(sig, "Unexpected relation",
             ErrType("Space relation found =",
                 info->thisSpace, SRELTYPES));
      break;
  }
  return (retval);
}


// ***********************************************************************

Space Space::operator[](int n)
{
  if ((n < 0) || (n >= info->cpSize)) {
    errh.ErrorExit("Space Space::operator[](int)",
                   "n is out of range",
           ErrVal("n = ", n), *this);
  }
  if (info->cpSize == 1) {
    return *this;
  } else {
    return(info->cpSpaces[n]);
  }
}

// ***********************************************************************
// ***********************************************************************
//
// VSpace Class
//
// ***********************************************************************
// ***********************************************************************
//
// Private/protected member functions
//
// ***********************************************************************
//
// This private method is responsible for system-built dual and tangent
// spaces.

VSpace::VSpace(SRel s, Space& ins)
{
  Space temp;
  GeObType nativetype;
  SRel access;
  Boolean useDual;
  char* tag;

// Catch the case where we are trying to build the dual of the Real
// Space, which is self-dual!  All other spaces Rn are distinct
// fom their dual spaces.

  if (s == LIN_DUAL) {
    temp = ins.GetSpace(LINEARIZATION);
    if (temp == Reals) {
      temp.SetDual(temp);
      temp.SetSpaceSet(LIN_DUAL, temp);
      temp.MergeSets(ins);
      return;
    }
  }

// Handle the regular cases.  Allocate a new SpaceInfo:

  info = new SpaceInfo;
  type = VEC_SPACE;
  holds = VEC_SPACE;

// Some stuff depends on what kind of space we are building:

  if (s == LIN_DUAL) {
    useDual = TRUE;
    access = LINEARIZATION;
  } else if (s == TANG_DUAL) {
    useDual = TRUE;
    access = TANGENT;
  } else if (s == TANGENT) {
    useDual = FALSE;
    access = AFFINE;
  } else {
    errh.ErrorExit("VSpace::VSpace(SRel, Space&)", "Unexpected relation",
           ErrType("Space relation found =", s, SRELTYPES));
  }

  temp = ins.GetSpace(access);

  if (useDual) {
    tag = this->BuildTagName("Dual to ", temp);
    nativetype = VECTOR;
    info->dual = temp;
  } else {
    tag = this->BuildTagName("Tangent to ", temp);
    nativetype = AFF_VECTOR;
  }

  this->CopyDebugName(tag);
  this->DestroyTagName(tag);

// Start to fill in all the information:

  info->type = VEC_SPACE;
  info->dim = temp.Dim();
  info->matsize = info->dim;
  info->native = nativetype;
  info->thisSpace = s;
  info->spaceSet[(int)s] = *this;
  info->stdBasis = Basis(VBASIS, *this, info->dim);
  info->isEuclidean = temp.IsEuclidean();

// If the given space is a CPSpace of spaces, this space should be
// a CPSpace of spaces, either the tangents or duals:

  info->cpSize = temp.CPSpaceSize();
  int dim_sum = 0;

  if (info->cpSize > 1) {
    info->cpSpaces = SpaceList(info->cpSize);
    info->cpSlots = IntList(info->cpSize);
    for (int j = 0; j < info->cpSize; j++) {
      if (useDual) {
        (info->cpSpaces)[j] = (temp[j]).Dual();
      } else {
        (info->cpSpaces)[j] = (temp[j]).GetSpace(TANGENT);
      }
      (info->cpSlots)[j] = dim_sum;
      dim_sum += (info->cpSpaces)[j].Dim();
    }
  } else {
    info->cpSlots = IntList(1);
  }

// Broadcast the space:

  this->Broadcast(s, ins);

// Build standard subsets:

  info->fullset = VSubSet(*this);
  info->fullproj = PSubSet(*this);

// Clean-up work (announcing dual or installing standard maps):

  if (useDual) {
    temp.SetDual(*this);
  } else {
    if (this->HasSpace(LINEARIZATION)) {
      this->LinkLT();
    }
  }
}

// ***********************************************************************
//
// Public member functions
//
// ***********************************************************************
//
// Used to cast a general space down to a vector space.
//

VSpace::VSpace(Space& s) : (s)
{
  type = VEC_SPACE;
  if ((holds != VEC_SPACE) && (holds != NULL_SPACE))  {
      errh.ErrorExit("VSpace::VSpace(Space&)",
                     "Attempted to cast a non-vector space to a vector space",
             s);
  }
}

// ***********************************************************************
//
// Memberwise assignment:
//

VSpace& VSpace::operator=(VSpace &s)
{
  holds = s.holds;
  info = s.info;
  return (*this);
}

// ***********************************************************************
//
// Output for debugging
//

void VSpace::debug_out(ostream &c, int indent)
{
  static char* name = "VSpace Object\n";
  this->data_out(c, indent, name);
  return;
}

// ***********************************************************************
//
// Stitches the given space into this space (which is a linearization) using
// the user specified map:
//

Space VSpace::SetSpace(Space& ins, Map& m)
{
  static char* sig = "Space VSpace::SetSpace(Space&, Map&)";

// Start by enforcing restrictions on the spaces.

  if (this->ThisSpaceIs() != LINEARIZATION) {
      errh.ErrorExit(sig,
                     "Only linearization spaces can be stitched manually",
             *this);
  }

  SRel rel = ins.ThisSpaceIs();

  if ((rel != AFFINE) && (rel != PROJECT_COMP)) {
    errh.ErrorExit(sig, "Incorrect type of space specified", ins);
  }

// The specified space must be unattached to another linear space, and
// we must be unattached:

  if (ins.HasSpace(LINEARIZATION)) {
    errh.ErrorExit(sig, "Specified space already linked", ins);
  }

  if (this->HasSpace(rel)) {
    errh.ErrorExit(sig, "This space already linked", *this);
  }

// Dimensionality must match:

  if (this->Dim() != ins.Dim() + 1) {
    errh.ErrorExit(sig, "Dimensions do not match", ins, *this);
  }

// If the other link of one is filled, the same link must be unfilled
// on the other.  We are taking advantage here of the fact that this
// space has not yet been registered with others in the set:

  if ((ins.HasSpace(AFFINE)) && (this->HasSpace(AFFINE))) {
    errh.ErrorExit(sig, "Linked to mismatched affine spaces", ins, *this);
  }
  if ((ins.HasSpace(PROJECT_COMP)) && (this->HasSpace(PROJECT_COMP))) {
    errh.ErrorExit(sig, "Linked to mismatched proj. completions", ins, *this);
  }

// Now check that the map we have been given fits the bill for stitching
// two spaces together.  If the given space is affine, the map must be
// affine; if projective, the map must be projective:

  if ((rel == AFFINE) && (m.Holds() != AFF_MAP)) {
    errh.ErrorExit(sig, "Map must be affine", ins, *this, m);
  }

  if ((rel == PROJECT_COMP) && (m.Holds() != PROJ_MAP)) {
    errh.ErrorExit(sig, "Map must be projective", ins, *this, m);
  }

// The range and domain must match the spaces we have:

  SubSet Rsub = m.Range();
  SubSet Dsub = m.Domain();
  Space Remb = Rsub.EmbeddingSpace();
  Space Demb = Dsub.EmbeddingSpace();

  if ((Remb != ins) || (!Rsub.IsFullSpace())) {
    errh.ErrorExit(sig, "Map Range/Dimension mismatch", ins, *this, m);
  }
  if ((rel == AFFINE) &&
      ((Demb != *this) || (Dsub.Dim() != ins.Dim()))) {
    errh.ErrorExit(sig, "Map Domain/Dimension mismatch", ins, *this, m);
  }
  if ((rel == PROJECT_COMP) &&
      ((Demb != *this) || (Dsub.Dim() != ins.Dim() + 1))) {
    errh.ErrorExit(sig, "Map Domain/Dimension mismatch", ins, *this, m);
  }

// Install the standard linkage maps.  Do the set merging only after
// we have checked the state of the original sets for choosing:

  if (rel == AFFINE) {
    if (this->HasSpace(PROJECT_COMP)) {
      this->MergeSets(ins);
      this->LinkLAHaveLP(m);
    } else if (ins.HasSpace(PROJECT_COMP)) {
      this->MergeSets(ins);
      this->LinkLAHaveAP(m);
    } else {
      this->MergeSets(ins);
      this->LinkLA(m);
    }
  } else {  // rel == PC
    if (this->HasSpace(AFFINE)) {
      this->MergeSets(ins);
      this->LinkLPHaveLA(m);
    } else if (ins.HasSpace(AFFINE)) {
      this->MergeSets(ins);
      this->LinkLPHaveAP(m);
    } else {
      this->MergeSets(ins);
      this->LinkLP(m);
    }
  }
  return (*this);
}

// ***********************************************************************
//
// Construct a plain vanilla Vector space, with no attachments to
// other spaces:
//

VSpace::VSpace(char* namein, int n, Boolean eflag)
{
  static char* sig = "VSpace::VSpace(char*, int, Boolean)";

  if (n < 1) {
    errh.ErrorExit(sig, "Dimension is out of range", ErrVal("Dim = ", n));
  }

// Create a new SpaceInfo:

  info = new SpaceInfo;
  type = VEC_SPACE;
  holds = VEC_SPACE;

// Fill in all the pertinent information:

  info->type = type;
  info->dim = n;
  info->matsize = info->dim;
  info->native = VECTOR;
  info->cpSize = 1;
  info->cpSlots = IntList(1);
  info->thisSpace = LINEARIZATION;
  info->stdBasis = Basis(VBASIS, *this, n);
  info->isEuclidean = eflag;
  this->CopyDebugName(namein);

// Fill in this slot of the space set:

  info->spaceSet[LINEARIZATION] = *this;

// Build standard subsets:

  info->fullset = VSubSet(*this);
  info->fullproj = PSubSet(*this);
}

// ***********************************************************************
//
// Build a vector space that is attached to specified affine and
// projective spaces.  Stitch them together with standard maps:
//

VSpace::VSpace(char* namein, Boolean eflag, ASpace& a, PSpace& p)
{
  static char* sig = "VSpace::VSpace(char*, Boolean, ASpace&, PSpace&)";

  type = VEC_SPACE;

// Make sure that the various spaces specified have the correct dimension:

  if (a.Dim() != p.Dim()) {
    errh.ErrorExit(sig, "Mismatch of dimensions", a, p);
  }

// The specified spaces must be unattached to linearizations:

  if (a.HasSpace(LINEARIZATION) || p.HasSpace(LINEARIZATION)) {
    errh.ErrorExit(sig, "One or more specified spaces already linked", a, p);
  }

// Start by constructing an unattached vector space:

  *this = VSpace(namein, a.Dim() + 1, eflag);

// Figure out if the specified spaces are unlinked or linked to each other,
// then link them if necessary:

  Boolean linkerr;
  Boolean linked = this->HaveLink(a, AFFINE, p, PROJECT_COMP, &linkerr);

  if (linkerr) {
    errh.ErrorExit(sig, "Spaces must be unlinked or linked to each other",
           a, p);
  }

  if (!linked) {
    p.MergeSets(a);
    p.LinkAP(AffineMap(a, p));
  }

// Broadcast the space:

  this->Broadcast(LINEARIZATION, a);

// Install the standard maps:

  LinkLAHaveAP(AffineMap(*this, a));
}

// ***********************************************************************
//
// Build a vector space that is attached to the specified space, which
// must be either affine or projective.  Stitch them together with
// standard maps:
//

VSpace::VSpace(char* namein, Boolean eflag, Space& ins)
{
  static char* sig = "VSpace::VSpace(char*, Boolean, Space&)";

  type = VEC_SPACE;

// The specified space must be the right kind:

  SRel rel = ins.ThisSpaceIs();

  if ((rel != AFFINE) && (rel != PROJECT_COMP)) {
    errh.ErrorExit(sig, "Incorrect type of space specified", ins);
  }

// The specified space must be unattached to another linearization:

  if (ins.HasSpace(LINEARIZATION)) {
    errh.ErrorExit(sig, "Specified space already linked", ins);
  }

// Start by constructing an unattached vector space:

  *this = VSpace(namein, ins.Dim() + 1, eflag);

// Broadcast the space:

  this->Broadcast(LINEARIZATION, ins);

// Install the standard maps:

  if (rel == AFFINE) {
    if (ins.HasSpace(PROJECT_COMP)) {
      LinkLAHaveAP(AffineMap(*this, ins));
    } else {
      LinkLA(AffineMap(*this, ins));
    }
  } else { // rel = PROJECT_COMP
    if (ins.HasSpace(AFFINE)) {
      LinkLAHaveAP(AffineMap(*this, ins.GetSpace(AFFINE)));
    } else {
      LinkLP(ProjectiveMap(*this, ins));
    }
  }
}

// ***********************************************************************
//
// This constructor is used to build Cartesian Product spaces from
// a list of existing spaces.  The space that results is a bare space;
// no other spaces are linked to it.
//

VSpace::VSpace(char* namein, SpaceList& t, Boolean eflag)
{
  static char* sig = "VSpace::VSpace(char*, SpaceList&, Boolean)";
  int i;

// Error if no spaces in the list, or if there is only one space
// in the list (cannot create aliases for atomic spaces):

  if (t.Length() <= 1) {
    errh.ErrorExit(sig, "Not enough spaces in space list", t);
  }

// Crete a new SpaceInfo:

  info = new SpaceInfo;
  type = VEC_SPACE;
  holds = VEC_SPACE;

// All the spaces in the list must be vector spaces.  The dimension
// of this space is the sum of the dimensions of the spaces in the
// list.  The cardinality of this space is the sum of the
// cardinalities of the spaces in the list (spaces in the list
// may themselves be CPSpaces).  A problem with this system (which
// lacks global knowledge of CP spaces in the system is that the
// user (and the system) can currently create aliases of certain
// special system-built CP spaces (the tangent and dual spaces
// of CPSpaces).  These aliases will always sit as linearization
// spaces in a six-space set.

  info->dim = 0;
  info->cpSize = 0;

  for (i = 0; i < t.Length(); i++) {
    Space temp = t[i];
    if (temp.Holds() != VEC_SPACE) {
      errh.ErrorExit(sig, "Non-vector space in space list", t);
    }
    info->dim += temp.Dim();
    info->cpSize += temp.CPSpaceSize();
  }

// Build the necessary lists and go through the space list once more,
// filling in as we go:

  info->cpSpaces = SpaceList(info->cpSize);
  info->cpSlots = IntList(info->cpSize);
  int slot = 0;
  int dim_sum = 0;

  for (i = 0; i < t.Length(); i++) {
    Space temp = t[i];
    for (int j = 0; j < temp.CPSpaceSize(); j++) {
      (info->cpSpaces)[slot] = temp[j];
      (info->cpSlots)[slot++] = dim_sum;
      dim_sum += (temp[j]).Dim();
    }
  }

// This is a bare vector space:

  info->type = type;
  info->matsize = info->dim;
  info->native = VECTOR;
  info->thisSpace = LINEARIZATION;
  info->stdBasis = Basis(VBASIS, *this, info->dim);
  info->spaceSet[LINEARIZATION] = *this;
  info->isEuclidean = eflag;
  this->CopyDebugName(namein);

// Build standard subsets:

  info->fullset = VSubSet(*this);
  info->fullproj = PSubSet(*this);

}

// ***********************************************************************
// ***********************************************************************
//
// ASpace Class
//
// ***********************************************************************
// ***********************************************************************
// ***********************************************************************
//
// Private/protected member functions
//
// ***********************************************************************
//
// Private method used by ASpaces to build their hoist matrices:
//

Matrix ASpace::BuildHoist(int n)
{
// The hoist matrix converts the matrix rep. of vectors in the tangent
// space to the matrix rep. in the affine space.  It is just an
// identity matrix with an extra column of zeros.

  Matrix retval = ZeroMatrix(n, n + 1);
  for (int i = 0; i < n; i++) {
    retval[i][i] = 1.0;
  }
  return (retval);
}

// ***********************************************************************
//
// Public member functions
//
// ***********************************************************************
//
// Downcasts a general space to an affine space (if it is one):

ASpace::ASpace(Space& s) : (s)
{
  type = AFF_SPACE;
  if ((holds != AFF_SPACE) && (holds != NULL_SPACE)) {
      errh.ErrorExit("ASpace::ASpace(Space&)",
                     "Attempted to cast a non-affine space to an affine space",
             s);
  }
}

// ***********************************************************************
//
// Assignment
//

ASpace& ASpace::operator=(ASpace &s)
{
  holds = s.holds;
  info = s.info;
  return (*this);
}

// ***********************************************************************
//
// Output for debugging:
//

void ASpace::debug_out(ostream &c, int indent)
{
  static char *name = "ASpace Object\n";
  this->data_out(c, indent, name);
  return;
}

// ***********************************************************************
//
// Stitches the given space to this space (which is affine) using
// the user specified map:
//

Space ASpace::SetSpace(Space& ins, Map& m)
{
  static char* sig = "Space ASpace::SetSpace(Space&, Map&)";

// Start by enforcing restrictions on the spaces.

  SRel rel = ins.ThisSpaceIs();

  if ((rel != LINEARIZATION) && (rel != PROJECT_COMP)) {
    errh.ErrorExit(sig, "Incorrect type of space specified", ins);
  }

// The specified space must be unattached to another affine space, and
// we must be unattached:

  if (ins.HasSpace(AFFINE)) {
    errh.ErrorExit(sig, "Specified space already linked", ins);
  }

  if (this->HasSpace(rel)) {
    errh.ErrorExit(sig, "This space already linked", *this);
  }

// Dimensionality must match:

  int targetdim = ((rel == LINEARIZATION) ? ins.Dim() - 1 : ins.Dim());
  if (this->Dim() != targetdim) {
    errh.ErrorExit(sig, "Dimensions do not match", ins, *this);
  }

// If the other link of one is filled, the same link must be unfilled
// on the other:

  if ((ins.HasSpace(LINEARIZATION)) && (this->HasSpace(LINEARIZATION))) {
    errh.ErrorExit(sig, "Linked to mismatched linearizations", ins, *this);
  }
  if ((ins.HasSpace(PROJECT_COMP)) && (this->HasSpace(PROJECT_COMP))) {
    errh.ErrorExit(sig, "Linked to mismatched proj. completions", ins, *this);
  }

// Now check that the map we have been given fits the bill for stitching
// two spaces together.  The map must be affine:

  if (m.Holds() != AFF_MAP) {
    errh.ErrorExit(sig, "Map must be affine", ins, *this, m);
  }

// The range and domain must match the spaces we have:

  SubSet Rsub = m.Range();
  SubSet Dsub = m.Domain();
  Space Remb = Rsub.EmbeddingSpace();
  Space Demb = Dsub.EmbeddingSpace();

  if ((Remb != ins) || (Rsub.Dim() != this->Dim())) {
    errh.ErrorExit(sig, "Map Range/Dimension mismatch", ins, *this, m);
  }
  if ((Demb != *this) || (!Dsub.IsFullSpace())) {
    errh.ErrorExit(sig, "Map Domain/Dimension mismatch", ins, *this, m);
  }

// Install the standard linkage maps.  Do the set merging only after
// we have checked the state of the original sets for choosing:

  if (rel == LINEARIZATION) {
    if (ins.HasSpace(PROJECT_COMP)) {
      this->MergeSets(ins);
      this->LinkLAHaveLP(m.Inv());
    } else if (this->HasSpace(PROJECT_COMP)) {
      this->MergeSets(ins);
      this->LinkLAHaveAP(m.Inv());
    } else {
      this->MergeSets(ins);
      this->LinkLA(m.Inv());
    }
  } else {  // rel == PC
    if (ins.HasSpace(LINEARIZATION)) {
      this->MergeSets(ins);
      this->LinkAPHaveLP(m);
    } else if (this->HasSpace(LINEARIZATION)) {
      this->MergeSets(ins);
      this->LinkAPHaveLA(m);
    } else {
      this->MergeSets(ins);
      this->LinkAP(m);
    }
  }
  return (*this);
}

// ***********************************************************************
//
// Construct a plain vanilla Affine space, with no attachments to
// other spaces:
//

ASpace::ASpace(char* namein, int n, Boolean eflag)
{
  static char* sig = "ASpace::ASpace(char*, int, Boolean)";

  if (n < 1) {
    errh.ErrorExit(sig, "Dimension is out of range", ErrVal("Dim = ", n));
  }

// Create a new SpaceInfo:

  info = new SpaceInfo;
  type = AFF_SPACE;
  holds = AFF_SPACE;

// Start to fill in the information:

  info->type = type;
  info->dim = n;
  info->matsize = info->dim + 1;
  info->native = AFF_POINT;
  info->thisSpace = AFFINE;
  info->stdBasis = Basis(FRAME, *this, n + 1);
  info->cpSize = 1;
  info->cpSlots = IntList(1);
  info->spaceSet[AFFINE] = *this;
  info->isEuclidean = eflag;
  this->CopyDebugName(namein);

// Need to construct the matrix that hoists tangent vectors into the
// internal affine space representation:

  info->hoistTangent = this->BuildHoist(info->dim);

// Build standard subsets:

  info->fullset = ASubSet(*this);
  info->stdAffineSet = info->fullset;
}

// ***********************************************************************
//
// Build an affine space that is attached to specified vector and
// projective spaces.  Stitch them together with standard maps:
//

ASpace::ASpace(char* namein, Boolean eflag, VSpace& v, PSpace& p)
{
  static char* sig = "ASpace::ASpace(char*, Boolean, VSpace&, PSpace&)";

  type = AFF_SPACE;

// Make sure that the various spaces specified have the correct dimension:

  if (v.Dim() != (p.Dim() + 1)) {
    errh.ErrorExit(sig, "Mismatch of dimensions", v, p);
  }

// The vector space must be a Linearization, not a dual or tangent space:

  if (v.ThisSpaceIs() != LINEARIZATION) {
    errh.ErrorExit(sig, "Vector space not a linearization", v);
  }

// The specified spaces must not be attached to an affine already:

  if (v.HasSpace(AFFINE) || p.HasSpace(AFFINE)) {
    errh.ErrorExit(sig, "One or more specified spaces already linked", v, p);
  }

// Start by constructing an unattached affine space:

  *this = ASpace(namein, p.Dim(), eflag);

// Figure out if the specified spaces are unlinked or linked to each other,
// then link them if necessary:

  Boolean linkerr;
  Boolean linked = this->HaveLink(v, LINEARIZATION, p, PROJECT_COMP, &linkerr);

  if (linkerr) {
    errh.ErrorExit(sig, "Spaces must be unlinked or linked to each other",
           v, p);
  }

  if (!linked) {
    p.MergeSets(v);
    p.LinkLP(ProjectiveMap(v, p));
  }

// Broadcast this space:

  this->Broadcast(AFFINE, v);

// Install the standard maps:

  LinkLAHaveLP(AffineMap(v, *this));
}

// ***********************************************************************
//
// Build an affine space that is attached to the specified space, which
// must be either vector or projective.  Stitch them together with
// standard maps:
//

ASpace::ASpace(char* namein, Boolean eflag, Space& ins)
{
  static char* sig = "ASpace::ASpace(char*, Boolean, Space&)";

  type = AFF_SPACE;

// The specified space must be the right kind:

  SRel rel = ins.ThisSpaceIs();

  if ((rel != LINEARIZATION) && (rel != PROJECT_COMP)) {
    errh.ErrorExit(sig, "Incorrect type of space specified", ins);
  }

// The specified space must be unattached to another affine space.

  if (ins.HasSpace(AFFINE)) {
    errh.ErrorExit(sig, "Specified space already linked", ins);
  }

// Start by constructing an unattached affine space:

  int targetdim = ((rel == LINEARIZATION) ? ins.Dim() - 1 : ins.Dim());
  *this = ASpace(namein, targetdim, eflag);

// Broadcast the space:

  this->Broadcast(AFFINE, ins);

// Install the standard maps:

  if (rel == LINEARIZATION) {
    if (ins.HasSpace(PROJECT_COMP)) {
      LinkLAHaveLP(AffineMap(ins, *this));
    } else {
      LinkLA(AffineMap(ins, *this));
    }
  } else { // rel = PROJECT_COMP
    if (ins.HasSpace(LINEARIZATION)) {
      LinkLAHaveLP(AffineMap(ins.GetSpace(LINEARIZATION), *this));
    } else {
      LinkAP(AffineMap(*this, ins));
    }
  }
}

// ***********************************************************************
//
// This constructor is used to build Cartesian Product spaces from
// a list of existing spaces.  The space that results is a bare space;
// no other spaces are linked to it.
//

ASpace::ASpace(char *namein, SpaceList& t, Boolean eflag)
{
  int i;

  static char* sig = "ASpace::ASpace(char*, SpaceList&, Boolean)";

// Error if no spaces in the list, or if one space in list (cannot
// create aliases for atomic spaces):

  if (t.Length() <= 1) {
    errh.ErrorExit(sig, "Not enough spaces in space list", t);
  }

// Create a new SpaceInfo:

  info = new SpaceInfo;
  type = AFF_SPACE;
  holds = AFF_SPACE;

// All the spaces in the list must be affine spaces.  The dimension
// of this space is the sum of the dimensions of the spaces in the
// list.  The cardinality of this space is the sum of the
// cardinalities of the spaces in the list (spaces in the list
// may themselves be CPSpaces).

  info->dim = 0;
  info->cpSize = 0;
  for (i = 0; i < t.Length(); i++) {
    Space hold = t[i];
    if (hold.Holds() != AFF_SPACE) {
      errh.ErrorExit(sig, "Non-affine space in space list", t);
    }
    info->dim += hold.Dim();
    info->cpSize += hold.CPSpaceSize();
  }

// Build the necessary lists and go through the space list once more,
// filling in as we go:

  info->cpSpaces = SpaceList(info->cpSize);
  info->cpSlots = IntList(info->cpSize);
  int slot = 0;
  int dim_sum = 0;

  for (i = 0; i < t.Length(); i++) {
    Space hold = t[i];
    for (int j = 0; j < hold.CPSpaceSize(); j++) {
      (info->cpSpaces)[slot] = hold[j];
      (info->cpSlots)[slot++] = dim_sum;
      dim_sum += (hold[j]).Dim();
    }
  }

// Misc. info:

  info->matsize = info->dim + 1;
  info->native = AFF_POINT;
  info->type = type;
  info->thisSpace = AFFINE;
  info->stdBasis = Basis(FRAME, *this, info->dim + 1);
  info->spaceSet[AFFINE] = *this;
  info->isEuclidean = eflag;
  this->CopyDebugName(namein);

// Need to construct the matrix that hoists tangent vectors into the
// internal affine space representation:

  info->hoistTangent = this->BuildHoist(info->dim);

// Build standard subsets:

  info->fullset = ASubSet(*this);
  info->stdAffineSet = info->fullset;
}

// ***********************************************************************
// ***********************************************************************
//
// PSpace Class
//
// ***********************************************************************

// ***********************************************************************
//
// Public member functions
//
// ***********************************************************************
//
// Used for downcasting:
//

PSpace::PSpace(Space& s) : (s)
{
  type = PROJ_SPACE;
  if ((holds != PROJ_SPACE) && (holds != NULL_SPACE)) {
      errh.ErrorExit("PSpace::PSpace(Space&)",
                     "Attempted to cast a non-proj. space to a proj. space",
             s);
  }
}

// ***********************************************************************
//
// Assignment:
//

PSpace& PSpace::operator=(PSpace &s)
{
  holds = s.holds;
  info = s.info;
  return (*this);
}

// ***********************************************************************
//
// Output for debugging:
//

void PSpace::debug_out(ostream &c, int indent)
{
  static char* name = "PSpace Object\n";
  this->data_out(c, indent, name);
  return;
}

// ***********************************************************************
//
// Stitches the given space into this space using the user specified map:
//

Space PSpace::SetSpace(Space& ins, Map& m)
{
  static char* sig = "void PSpace::SetSpace(Space&, Map&)";

// Start by enforcing restrictions on the spaces.
// The specified space must be the right kind:

  SRel rel = ins.ThisSpaceIs();

  if ((rel != AFFINE) && (rel != LINEARIZATION)) {
    errh.ErrorExit(sig, "Incorrect type of space specified", ins);
  }

// The specified space must be unattached to another PC space, and
// we must be unattached:

  if (ins.HasSpace(PROJECT_COMP)) {
    errh.ErrorExit(sig, "Specified space already linked", ins);
  }

  if (this->HasSpace(rel)) {
    errh.ErrorExit(sig, "This space already linked", *this);
  }

// Dimensionality must match:

  int targetdim = ((rel == LINEARIZATION) ? ins.Dim() - 1 : ins.Dim());
  if (this->Dim() != targetdim) {
    errh.ErrorExit(sig, "Dimensions do not match", ins, *this);
  }

// If the other link of one is filled, the same link must be unfilled
// on the other:

  if ((ins.HasSpace(AFFINE)) && (this->HasSpace(AFFINE))) {
    errh.ErrorExit(sig, "Linked to mismatched affine spaces", ins, *this);
  }
  if ((ins.HasSpace(LINEARIZATION)) && (this->HasSpace(LINEARIZATION))) {
    errh.ErrorExit(sig, "Linked to mismatched linearizations", ins, *this);
  }

// Now check that the map we have been given fits the bill for stitching
// two spaces together.  If the given space is affine, the map must be
// affine; if linear, the map must be projective:

  if ((rel == AFFINE) && (m.Holds() != AFF_MAP)) {
    errh.ErrorExit(sig, "Map must be affine", ins, *this, m);
  }

  if ((rel == LINEARIZATION) && (m.Holds() != PROJ_MAP)) {
    errh.ErrorExit(sig, "Map must be projective", ins, *this, m);
  }

// The range and domain must match the spaces we have:

  SubSet Rsub = m.Range();
  SubSet Dsub = m.Domain();
  Space Remb = Rsub.EmbeddingSpace();
  Space Demb = Dsub.EmbeddingSpace();

  if ((Remb != ins) || (!Rsub.IsFullSpace())) {
    errh.ErrorExit(sig, "Map Range/Dimension mismatch", ins, *this, m);
  }
  if ((rel == AFFINE) &&
      ((Demb != *this) || (Dsub.Dim() != ins.Dim()))) {
    errh.ErrorExit(sig, "Map Domain/Dimension mismatch", ins, *this, m);
  }
  if ((rel == LINEARIZATION) &&
      ((Demb != *this) || (Dsub.Dim() != ins.Dim() - 1))) {
    errh.ErrorExit(sig, "Map Domain/Dimension mismatch", ins, *this, m);
  }

// Install the standard linkage maps.  Do the set merging only after
// we have checked the state of the original sets for choosing:

  if (rel == AFFINE) {
    if (this->HasSpace(LINEARIZATION)) {
      this->MergeSets(ins);
      LinkAPHaveLP(m.Inv());
    } else if (ins.HasSpace(LINEARIZATION)) {
      this->MergeSets(ins);
      LinkAPHaveLA(m.Inv());
    } else {
      this->MergeSets(ins);
      LinkAP(m.Inv());
    }
  } else {  // rel == LINEARIZATION
    if (this->HasSpace(AFFINE)) {
      this->MergeSets(ins);
      LinkLPHaveAP(m.Inv());
    } else if (ins.HasSpace(AFFINE)) {
      this->MergeSets(ins);
      LinkLPHaveLA(m.Inv());
    } else {
      this->MergeSets(ins);
      LinkLP(m.Inv());
    }
  }
  return (*this);
}

// ***********************************************************************
//
// Construct a plain vanilla projective space, with no attachments to
// other spaces:
//

PSpace::PSpace(char *namein, int n)
{
  static char* sig = "PSpace::PSpace(char*, int)";

  if (n < 1) {
    errh.ErrorExit(sig, "Dimension is out of range", ErrVal("Dim = ", n));
  }

// Create a new SpaceInfo:

  info = new SpaceInfo;
  type = PROJ_SPACE;
  holds = PROJ_SPACE;

// Fill in all the pertinent information:

  info->type = type;
  info->dim = n;
  info->matsize = info->dim + 1;
  info->native = PROJ_POINT;
  info->thisSpace = PROJECT_COMP;
  info->stdBasis = Basis(HFRAME, *this, n + 1);
  info->spaceSet[PROJECT_COMP] = *this;
  info->isEuclidean = FALSE;
  info->cpSize = 1;
  info->cpSlots = IntList(1);
  this->CopyDebugName(namein);

// Build standard subsets:

  info->fullset = PSubSet(*this);
  info->fullproj = info->fullset;
}

// ***********************************************************************
//
// Build a projective space that is attached to specified vector and
// affine spaces.  Stitch them together with standard maps:
//

PSpace::PSpace(char* namein, VSpace& v, ASpace& a)
{
  static char* sig = "PSpace::PSpace(char*, VSpace&, ASpace&)";

  type = PROJ_SPACE;

// Make sure that the various spaces specified have the correct dimension:

  if (v.Dim() != (a.Dim() + 1)) {
    errh.ErrorExit(sig, "Mismatch of dimensions", v, a);
  }

// The vector space must be a Linearization, not a dual or tangent space:

  if (v.ThisSpaceIs() != LINEARIZATION) {
    errh.ErrorExit(sig, "Vector space not a linearization", v);
  }

// The specified spaces must be unattached to PCs:

  if (v.HasSpace(PROJECT_COMP) || a.HasSpace(PROJECT_COMP)) {
    errh.ErrorExit(sig, "One or more specified spaces already linked", v, a);
  }

// Start by constructing an unattached projective space:

  *this = PSpace(namein, a.Dim());

// Figure out if the specified spaces are unlinked or linked to each other,
// then link them if necessary:

  Boolean linkerr;
  Boolean linked = this->HaveLink(a, AFFINE, v, LINEARIZATION, &linkerr);

  if (linkerr) {
    errh.ErrorExit(sig, "Spaces must be unlinked or linked to each other",
           a, v);
  }

  if (!linked) {
    v.MergeSets(a);
    v.LinkLA(AffineMap(v, a));
  }

// Broadcast this space:

  this->Broadcast(PROJECT_COMP, a);

// Install the standard maps:

  LinkAPHaveLA(AffineMap(a, *this));
}

// ***********************************************************************
//
// Build a projective space that is attached to the specified space, which
// must be either vector or affine.  Stitch them together with
// standard maps:
//

PSpace::PSpace(char* namein, Space& ins)
{
  static char* sig = "PSpace::PSpace(char*, Space&)";

  type = PROJ_SPACE;

// The specified space must be the right kind:

  SRel rel = ins.ThisSpaceIs();

  if ((rel != LINEARIZATION) && (rel != AFFINE)) {
    errh.ErrorExit(sig, "Incorrect type of space specified", ins);
  }

// The specified space must be unattached to another projective space.

  if (ins.HasSpace(PROJECT_COMP)) {
    errh.ErrorExit(sig, "Specified space already linked", ins);
  }

// Start by constructing an unattached affine space:

  int targetdim = ((rel == LINEARIZATION) ? ins.Dim() - 1 : ins.Dim());
  *this = PSpace(namein, targetdim);

// Broadcast the space:

  this->Broadcast(PROJECT_COMP, ins);

// Install the standard maps:

  if (rel == AFFINE) {
    if (ins.HasSpace(LINEARIZATION)) {
      LinkAPHaveLA(AffineMap(ins, *this));
    } else {
      LinkAP(AffineMap(ins, *this));
    }
  } else { // rel = LINEARIZATION
    if (ins.HasSpace(AFFINE)) {
      LinkAPHaveLA(AffineMap(ins.GetSpace(AFFINE), *this));
    } else {
      LinkLP(ProjectiveMap(ins, *this));
    }
  }
}


// ***********************************************************************
// ***********************************************************************
//
// Create the error information block:
//

static SpaceInfo errInfo("Uninitialized Space", 0);

// ***********************************************************************
//
// Create the vector space for Real numbers:
//

VSpace Reals("Real Numbers", 1, FALSE);

// ***********************************************************************